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1 INTRODUCTION 



ABSTRACT 

We study the evolution of phase-space density during the hierarchical structure formation 
of ACDM halos. We compute both a spherically-averaged surrogate for phase-space density 
(Q = p/cr 3 ) and the coarse-grained distribution function /(x, v) for dark matter particles that 
lie within ~ 2 virial radii of four Milky- Way-sized dark matter halos. The estimated /(x, v) 
spans over four decades at any radius. Dark matter particles that end up within two virial 
radii of a Milky- Way-sized DM halo at z — have an approximately Gaussian distribution 
in log(/) at early redshifts, but the distribution becomes increasingly skewed at lower red- 
shifts. The value / pC ak corresponding to the peak of the Gaussian decreases as the evolution 
progresses and is well described by / pea k(z) oc (1 + z) for z > 1. The highest values of 
/, (responsible for the skewness of the profile) are found at the centers of dark matter halos 
and subhalos, where / can be an order of magnitude higher than in the center of the main 
halo. We confirm that Q(r) can be described by a power-law with a slope of —1.8 ± 0.1 over 
2.5 orders of magnitude in radius and over a wide range of redshifts. This Q(r) profile likely 



reflects the distribution of entropy (K = a 2 j p 



2/3 
dm 



,1.2 



), which dark matter acquires as it 



is accreted onto a growing halo. The estimated /(x, v), on the other hand, exhibits a more 
complicated behavior. Although the median coarse-grained phase-space density profile F(r) 
can be approximated by a power-law, oc r - 16 ±o.i5^ m tne inner regions of halos (< 0.6 r V j r ), 
at larger radii the profile flattens significantly. This is because phase-space density averaged 
on small scales is sensitive to the high-/ material associated with surviving subhalos, as well 
as relatively unmixed material (probably in streams) resulting from disrupted subhalos, which 
contribute a sizable fraction of matter at large radii. 



Cosmological iV-body simulations of the formation of structure in 
the Universe allow us to reconstruct how gravitationally bound ob- 
jects like galaxies and clusters form and evolve. These simulations 
have shown that despite the seemingly complex hierarchical for- 
mation history of dark matter (DM) halos, the matter density dis- 
tributions of cosmolo gical DM halos are described by si mple 3- 
parameter profiles dNavarro eTal . 1996; Merritt et al. 2006). 



The study of the evolution of the phase-space distribution 
function (DF) of a collisionless dynamical system (composed of 
stars or DM or both) is fundamental to understanding its dynami- 
cal evolution, since a collisionless system is completely described 
by its phase-space density distribution (otherwise called the fine- 
grained DF f (x, v)). For an isolated collisionless system, f (x, v, t) 
is fully described by the Collisionless Boltzmann Equation (CBE) 
(sometimes called the Vlasov equation), which is a continuity equa- 
tion in phase-space. A consequence of the CBE is that the fine- 
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grained DF is always conserved. In addition, V(f)df, the volume 
of phase-space occupied by phase-space elements whose density 
lies between (f , f + df), is also conserved. However, the con- 
servation of the fine-grained f(x, v) and V(f)df is not a very 
useful property for understanding the evolution of real systems, 
since in practice it is only possible to compute an average of f 
over a finite volume of phase-space. This average is referred to as 
the coarse-grained DF /(x, v), and the associated volume density 
is V(f). The evo lution of the coarse-grained DF is governed by 
Mixing Theor ems dBinnev & Tremain el 1987l : lTremaine et alj 19861 ; 
Mathu rTl988l) . These theorems state that if the coarse-grained DF is 
bounded, then processes that operate during the relaxation of colli- 
sionless systems (e.g. phase mixing, chaotic mixing, and the mixing 
of energy and angular momentum that accompanies violent relax- 
ation) result in a decrease in the coarse-grained phase-space den- 
sity. This is because at any point in phase-spac e, both low and high 
phase-space density regions get highly mixed dBinnev & Tremaind 
1 19871 : 1 Tremaine et al.|[l986f) . In addition, there are simple relation- 
ships that exist between the initial V(f) and the coarse- grained 
V(f) at a later time: V(f) is always greater than V(f) dMathuj 
Il988h . 

Recently, iDehnerJ d2005l) proved a new Mixing Theorem that 
states that the excess-mass function D(f), which is the mass of ma- 
terial with coarse-grained phase-space density greater than a value 
/, always decreases due to mixing, for all values of /. He showed, 
using this theorem, that steeper cusps are less mixed than shallower 
ones, independent of the details of the DF or density profile, and 
this implies that a merger remnant cannot have a cusp that is either 
steeper or shallower than the steepest of its progenitors. 

iTavlor & Navarrol d200lh have introduced a surrogate quan- 
tity, Q(r) — p(r)/a(r) 3 , where p(r) is the configuration-space 
density averaged in spherical shells, and a(r) is the velocity dis- 
persion of DM particles averaged in a spherical shell centered at 
the radius r. Note that, although Q(r) has dimensions of phase- 
space density, it is not true coarse-grained phase-space density, be- 
cause it is constructed out of separately computed configuration- 
space density and velocity dispersion in arbitrarily chosen spher- 
ical shells. Consequently, it is not expe cted to obey any Mixing 
Theorems (Dehnen & McLaughlin 20051) an d is generally difficult 
to interpret in terms of phas e-space density. Neverthel ess, it is an 
intriguing quantity, because ITavlor & Nava rro ( 2001) show that 
Q(r) is well approximated by a single power-law, Q cx r~ 13 with 
(3 ~ 1.874 over more than 2.5 orders of magnitude in radius for 
CDM halos univ ersally, regardless of mass and backg round cos- 
mology (see also lRasia et alj|2004l ; lAscasibar et al. 2004). Recent 
work, however, indicates that Q(r) profiles are somewhat sensitive 
to the amount of subs t ructure and the slope of the power spectrum 
dWang & Whitel2008l:lKnollmarm et aT]|2008h . 

While power-law Q(r) profiles are by no me ans universal 
for se lf-gravitating systems in dynamical equilibrium (Barn es et al.l 
2006), power-law profiles with exactly the same slope were first 
shown to arise in simple self-similar spherical gravitational infall 



models dBertschinger| [l985) pointing to a possible un iversality in 
the mechanisms th at produces th em. Austin et"al] d2005h extended 
the early work of iBertschinged d 19851) using semi-analytical ex- 
tended secondary infall models to follow the evolution of collision- 
less spherical shells of matter that are initially set to be out of dy- 
namical equilibrium and are allowed to move only radially. They 
concluded that the power-law behavior of the final phase-space 
density profile is a robust feature of virialized halos, which have 
reached equilibrium via violent relaxation. They used a constrained 
Jeans equation analysis to show that this equation has different 
types of solutions and that it admits a unique periodic solution 
which gives a power-law density profile with slope /3 — 1.9444 
(comparable to the num erical value of /3 = 1.87 obtained by 
ITavlor & Navarrol d200lh ). 

A recent study by iHoffman et alj d2007h followed the evo- 
lution of Q(r) profiles in cosmological DM halos in constrained 
simulations designed to control the merging history of a given 
halo. These authors showed that, during relatively quiescent phases 
of halo evolution, the density profile closely follows that of a 
NFW halo, and Q(r) is always well represented by a power-law 
r~P with P — 1.9 ± 0.1. They showed that Q{r) deviates from 
power-law most strongly during major mergers but recovers the 
power-law form thereafter (this is consistent wit h our findings us - 
ing a series of controlled m e rger s imulations, IVass et al 2009). 
More recently, IWoitak et alj d2008l) demonstrated (using a vari- 
ant of the FiEstAScode) that the DF of ACDM halos of mass 
10 14 - 10 15 M Q can be separated into energy and angular mo- 
mentum components. They proposed a phenomenological model 
for spherical potentials and showed that their model DF was a good 
match to the N-body DFs and was also able to reproduce the power- 
law behavior of Q(r). Another stu dy of relevance is the work of 
IPeirani & de Freitas Pachecol d2007h . which defined a global value 
for the phase-space density, Q, which they found decreases rapidly 
with time. Despite the insights obtained in these studies, the origin 
for the universality of Q(r) is not yet understood, and the ques- 
tion of how this quantity relates to the true coarse-grained phase- 
space density has not yet been thoroughly investigated. This moti- 
vates our study which is focused on the evolution of both Q(r) and 
phase-space density evolution using a set of self-consistent cosmo- 
logical simulations of halo formation. 

In the last four years, three independent numerical codes 
to compute the true coarse-grained phase-sp ace density from 

dXTad et all 120041: 



A' -body simulations have been de veloped ( Arad et al.l 120041 
lAscasibar & Binnevll2005l : ISharma & Steinmetzll2006l) . In the Ap 



pendix we compare results obtained with t he codes "FiEstAS" 
dAscas ibar & Binnev 2005) and "EnBiD" dSharma & Stei nmetz 
2006) at a range of redshifts. However we only present results 
obtained with the "EnBiD" code using a specific kernel. Follow- 
ing the submission of thi s paper we became aware of the work of 
Macieiews ki et alJd2008l) . These authors carried out a similar com- 
parison of various coarse-grained phase space density estimators 
for N-body simulations cosmological halos at z = 0. 
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In an accompanying paper dVass et alj2009h . we present anal- 
ysis of phase-space and Q(r) evolution during major mergers be- 
tween CDM-like spherical halos. We show that major mergers 
(those which are the most violent and therefore likely to result in 
the greatest amount of mixing in phase-space) preserve spherically- 
averaged phase-space density profiles (both Q(r) and F(r)) out to 
the v irial radius r v j r of the final remnant. We confirm the predic- 
tion IBehnen 2005) that the phase-space density profiles of DM 
halos are extremely robust, and in particular, the prediction that the 
steepest central cusp always survives. 

To understand th e robustness of profiles outsid e the central 
cusp dVass et alj |2009). we drew on recent work by IValluri et al.l 
(200 7) that shows that matter in equally-spaced radial shells is re- 
distributed during the merger is such a way that only about 20% of 
the matter in the central cusp is ejected during the merger to radii 
extending out to about two scale radii. For all other radii, roughly 
15% by mass in each radial shell is redistributed uniformly with 
radius. Further, nearly 40% of the mass of the pr ogenitor halos lies 
beyond the formal virial radius of the remnant dKazantzidis et al.l 
20061; IValluri et alj|2007h . and this matter originates roughly uni- 
formly from each radial interval starting from about three scale 
ra dii from the center. While major mergers (such as those studied 
bv lDehned B: IVass et alJ l^OOSh are instructive, they represent 
the most extreme form of mixing, and are not the major mode of 
mass accretion in the Universe. In addition, since experiments with 
major mergers demonstrate that pre-existing power-law profiles are 
preserved, they shed little light on the formation of the profiles. Our 
primary goal is to gain a better understanding of the origin of the 
power-law phase-space density profiles seen in cosmological DM 
simulations. 

In this paper we investigate the evolution of the coarse-grained 
phase-space density in the formation and evolution of four Milky- 
Way-sized halos in a ACDM cosmology with cosmological param- 
eters: (fi m , Qa, h, as) — (0.3, 0.7, 0.7, 0.9). In § [2] we summa- 
rize the cosmological N-body simulations analyzed in this study as 
well as the numerical methods used to obtain coarse-grained phase- 
space densities. (The Appendix contains a detailed comparison of 
the two codes FiEstAS and EnBiD . We also present reasons for our 
choice of EnBiD , with a n = 10 smoothing kernel). In §[3] we de- 
scribe the properties of the spherically-averaged quantity Q(r) as 
well as the the coarse-grained DF for one Milky-Way-sized DM 
halo at z = 0, as well as the evolution of the phase-space DF of 
this halo from z — 9 to the present. In §[4] we also present results 
for three additional Milky-Way-sized halos. In §[5] we present an 
interpretation of the power-law profiles seen in Q(r) and discuss 
the implications of observed coarse grained distribution function in 
a cosmological context. § [6] summarizes the results of this paper 
and concludes. 



2 NUMERICAL METHODS 

The simulations an alyzed in this paper ar e des cribed in greater de- 
tail in the works of iKravtsov et al.l d2004h and lGnedin & Kravtsov 
(2006). The simulations were car ried out using the Ad aptive Re- 
finement Tree Af-body code (ART IKravtsov et al.|[l997h . The sim- 
ulation starts with a uniform 256 3 grid covering the entire computa- 
tional box. This grid defines the lowest (zeroth) level of resolution. 
Higher force resolution is achieved in the regions corresponding to 
collapsing structures by recursive refining of all such regions by 
using an adaptive refinement algorithm. Each cell can be refined 
or de-refined individually. The cells are refined if the particle mass 
contained within them exceeds a certain specified value. The grid is 
thus refined to follow the collapsing objects in a quasi-Lagrangian 
fashion. 

Three of the galactic DM halos were simulated in a comoving 
box of 25 h^ 1 Mpc (hereafter L25); they were selected to reside 
in a well-defined filament at z = 0. Two halos are neighbors, lo- 
cated 425 h~ kpc from each other. The third halo is isolated and is 
located two Mpc away from the pair. Hereafter, we refer to the iso- 
lated halo as Gi and the halos in the pair as G2 and G3. Although 
we analyzed all three halos at the same level of detail we present 
detailed results for Gi . The virial masses and virial radi i for the ha- 
los studied are given in Table 1 in lKravtsov et al.l d2004l) . The virial 
radius (and the corresponding virial mass) was chosen as the radius 
encompassing a mean density of ~ 200 times the mean density of 
the Universe. The masses of the DM halos are well within the range 
of poss ibilities allowed by models for the halo of the Milky Way 
galaxy dKlypin et alj2002f) . The simulations followed a Lagrangian 
region corresponding to the sphere of radius equal to two virial radii 
around each halo at z = 0. This region was re-sampled with the 
highest resolution particles of mass m p = 1.2 x 10 6 /i~ 1 Mq, cor- 
responding to 1024 3 particles in the box, at the initial redshift of the 
simulation (z\ = 50). The maximum of ten refinement levels was 
reached in the simulations corresponding to the peak formal spa- 
tial resolution of 152/x" 1 parsec. Each host halo is resolved with 
~ 10 6 particles within its virial radius at z = 0. 

The fourth halo was simulated in a comoving box of 20 h~ l 
Mpc box (hereafter L20), and it was used to follow the Lagrangian 
region of approximately five virial radii around the Milky Way- 
sized halo with high resolution. In the high-resolution region the 
mass of the DM particles is m p = 6.1 x 10 s Mq, correspond- 
ing to effective 1024 3 particles in the box, at the initial redshift 
of the simulation (zi = 70). As in the other simulation, this run 
starts with a uniform 256 3 grid covering the entire computational 
box. Higher force resolution is achieved in the regions correspond- 
ing to collapsing structures by recursive refining of all such re- 
gions using an adaptive refinement algorithm. Only regions con- 
taining highest resolution particles were adaptively refined. The 
maximum of nine refinement levels was reached in the simulation 
corresponding to the peak formal spatial resolution of 150/i _1 co- 
moving parsec. The Milky Way-sized host halo has the virial mass 
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of 1.4 x 10 12 /i _1 M Q (or 2.3 million particles within the virial 
radius) and virial radius of 230/i -1 kpc at z — 0. 

The peak spatial resolution of the simulations determines the 
minimum radius to which we can trust the density and velocity pro- 
files. In the following analysis, we only consider profiles at radii at 
least eight times larger than the peak resolution of the simulations 
(i.e., « 0.5 — 0.6/i _1 comoving kpc or r « 0.004 r vir ( z =o))- 



2.1 Phase-space density estimators 

In this paper, we will focus on the time evolution of the spher- 
ically averaged phase-space density Q(r) as well as the coarse- 
grained phase-space DF /(x, v) during the formation of the Milky- 
Way sized DM halos in ACDM cosmological simulations described 
above. 

In the last four years three independent numerical codes 
to compute the coarse-grained phase-space density from N- 
body simulations have been developed dArad et all 1 20041 ; 
lAscasibar & Binnevl l2005t Isharma & Steinmetzl I2006HH These 
techniques differ in the scheme used to tesselate 6-dimensional 
phase space as well as in the density estimator s they use. The first 
study of coarse grained DF jArad et alj|2004h showed that it has 
its highest values at the centers of DM halos and subhalos. They 
also showed that the volume density of phase-space V(f) within 
individual cosmological DM halos at z = has a power-law 
profile over nearly four orders of magnitude in /. The main 
criticism of this code is that it is not metric-free. The FiEstAS 
al gorithm of lAscasibar & Binnevl d2005h and the EnBid algorithm 
of ISharma & Steinmetzl d2006h are preferred due to their speed 
and the metric free nature. The FiEstAS algorithm is based on 
a repeated division of each dimension of phase-space into two 
region s that contain roughly e qual numbers of particles. En- 
BiD {sh arma & Steinmetzl 2006) closely follows the method used 
in FiEstAS , but it optimizes the number of divisions to be made in 
a particular dimension by using a minimum-entropy criterion based 
on the Shannon entropy that measures the phase-space density 
with m uch greater accuracy when / is high. ISharma & Steinmetzl 
(2006) showed that EnBiD with a kernel which includes 10 nearest 
neighbors (n — 10) about each point was able to recover analytic 
phase-space density profiles to nearly 3-4 decades higher values of 
/ than FiEstAS . 



1 An ev en more sophisticated num erical method has recently been pre- 
sented by Vogelsberger et al. (2008) for calculating the fine-grained phase- 
space structure of DM distributions derived from cosmological simulations. 
This code has the potential to identify fine-scale structure such as caustics in 
phase-space and the phase-space structure of tidal streams in the Milky- Way 
halo. Its ability to estimate the fine-grained f (x, v) makes it useful for un- 
derstanding non-equilibrium systems and for making accurate predictions 
for direct DM searches. 
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Figure 1, Q(r) and F(r) in 100 spherical radial bins about the most- 
bound-particle in halo G\ from the L25 simulation at z=0. Top panel: Q 
(dot-dashed curve); power-law fit Q ex. r~ 1-84±0.012 j s gj ven by upper 
thin solid line. F(r) from EnBiD (with n = 10 smoothing kernel (dot- 
ted curve)); power-law fit to F(r) oc r~ 1 -59±0.054 j s gj ven i ower thin 
solid line. Mean value of f (long-dashed curve) is noisier than F(r) due 
to subhalos. Residuals to the two power-law fits are shown in the bottom 
panel. 



For our halos at z — 0, the results we obtained with the dif- 
ferent codes were complet ely consistent with those obtained by 
ISharma & Steinm etz (2006). The deviations between the estimates 
obtained with the different codes and various parameters were sig- 
nificantly higher at higher redshift - where comparisons with an- 
alytic estimates are not available. We refer the reader to the Ap- 
pendix for a detailed comparison between estimated values of / 
using the two codes at various redshifts. In the rest of this paper we 
present results o btained with EnBiD (n = 10 kernel), the parame- 
ters preferred by Sh arma & Steinmetzl d2006l) . 

Figure Q] Top shows the spherically averaged quantity Q = 
p/o- 3 (dot-dashed curve); the mean of /(x,v) in spherical bins 
(dashed curve) and the median DF in spherical bins (dotted curve). 
As we will show in Figure [4] the large fluctuations in the mean 
value of / (dashed curve) beyond 0.1rvi r are due to the presence 
of substructure, which can have extremely high central values of 
/. The median value of / (hereafter represented by F) is much 
smoother, being less sensitive to the large range (nearly eight or- 
ders of magnitude) in / at each radius. In what follows, we will 
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use the median F, computed in concentric radial bins centered on 
the most bound particle in the main halo, since it is less sensitive to 
substructure. Q(r ) is well fitted by a power-law: Q <x j.- 1 - 84 ^ 012 
over the radial range [0.004-0.6] r/r v i r (upper thin solid line). The 
power-law fit F(r) oc r - 1 - 59 ± 054 over th e same radial range is 
shown by the lower thin solid line. The bottom panel shows the 
residuals of the fits F (log(F/F fit )) and Q (log(Q/Q flt )). This 
plot shows that while Q(r) is an extremely good power-law, the 
median F(r) is only approximately power-law over the same ra- 
dial range. Note that Q(r) is numerically larger than F(r) due to 
the fact that the former quantity does not properly account for the 
volume of the phase-space element. A simi lar result was obtain ed 
in a much higher resolution simulation by Sta del et all d2008h . a 
study which appeared as we were preparing this paper for publica- 
tion. 



3 EVOLUTION OF PHASE-SPACE DENSITY WITH 
REDSHIFT 

Hoff man et alj d2007t) studied the phase-space density profiles of a 
DM halo by tracking Q(r) for material within the formal virial ra- 
dius at each redshift. They found that the virialized material within 
this radius has an approximately power-law form with a constant 
power-law index of —1.9 ± 0.05 at all redshifts from z = 5 to the 
present. 

We follow a slightly different approach in this paper, since our 
objective is to understand the evolution of the true coarse-grained 
phase-space density distribution. We track the evolution of phase- 
space density, by tracing backwards in time, all the material that lies 
inside the virial radius at z = 0. This will allow us to understand 
how the initially high phase-space density material that lies out- 
side virialized systems at high redshift falls in and undergoes mix- 
ing and how the resultant mixing preserves the phase-space density 
profiles as a function of redshift. 

We identified Milky-Way sized DM halos at z — in our cos- 
mological simulations and identify all the particles that lie inside 
twice the virial radius at z = 0. Our ch oice of halo outer radius 
(2r v ir) is motivated by the recent work o f ICuesta et al] d2008h who 
showed that the virial radius is arbitrary and does not correspond 
to a physically meaningful outer boundary for a DM halo. A better 
outer radius is the so-called "static radius", defined as the radius 
at which the mean radial velocity of particles is zero. This radius 
is typically about 2r v i r for a galaxy-sized DM halo. After identi- 
fying all particles within 2r v i r at z = 0, we tracked them back to 
z = 9. Our simulations do proceed back in time to even higher 
redshifts, however we do not analyze them here because the mass 
resolution at higher redshifts does not allow us to resolve the evolu- 
tion of most of the objects beyond this epoch. There were 1.6x 10 6 
particles within 2 virial radii of the Gi halo, which is the halo that 
is most isolated at z = 0. Two other halos (G2 and G3) from the 
L25 simulation are closer to each other at 2 = (1.5 virial radii 
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Figure 2. F as a function of comoving radius r for all material that lies 
within 2r v ; r of the center of the Gi halo at z = 0. Each curve corresponds 
to a different redshift as indicated in the legends. 



apart). Therefore, we only tracked those particles which lay within 
one virial radius of the centers of these two halos at z = back 
in time to z = 9. The particles in the high resolution simulation 
L20 were selected and tracked in the same way as for the Gi halo, 
and at z = there were 3x 10 6 particles inside twice the virial ra- 
dius. The position and velocity data for all the particles identified 
as belonging to a given halo at z — were analysed. 

At all redshifts we compute the phase-space densities of parti- 
cles in physical coordinates rather than in co-moving coordinates. 
Distances are always in kpc and velocities are in kms~ , and these 
quantities are measured relative to the center of the most massive 
progenitor. In many of the figures that follow, the physical distances 
at all redshifts other than z — are scaled by the virial radius at 
z = 0; i.e. all distances are given in units of r/r v i r ( z= n). 

In this section, we will present the results obtained for the Gi 
halo in the L25 simulation. A comparison with results for the other 
three halos is made in §|4] 

Figure [2] plots the median value (F) of the phase space den- 
sity / as a function of radius (as defined above) at several different 
redshifts. At the highest redshift plotted (z — 9), the overall value 
of F is higher at large radii than it is at the center of the halo, 
and it only varies by a factor of few over the entire range of radii 
(F ~ 10 - 3OM kpc" 3 (kms -1 ) -3 ). This is because within the 
inner, virialized regions the coarse-grained phase-space density is 
lowered due to mixing, while it is still high in the outer regions 
where the large fraction of matter is not yet mixed or is located in 
small subhalos, that have not mixed to the same degree as the main 
host. As the evolution progresses, there is a decrease in the central 
value of F till z = 3. The median value of F at the center of the 
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halo drops from F ~ lOM0kpc _3 (kms _1 )~ 3 at z — 9 down to 
F « 2M Q kpc~ 3 (kms _1 )~ 3 at 2 = 0. The central value of F 
remains constant beyond z — 3, while there is a steady decrease in 
F at intermediate radii. 

Although it is instructive to plot the mean and/or median val- 
ues of F at each radius, much more information is contained in 
the full coarse-grained phase-space density /. The six-dimensional 
function /(x, v) can be most easily visualized using the volum e 
DF V(.f), which also obeys a Mixing Theorem dMathuj[l988h . 
lArad et all (2004) showed that at z = 0, V(f) for cosmological 
DM halos follows a power-law profile V(f) cx f~ a with power- 
law index a — 2.4 over four orders of magnitude in /. They argued 
that since DM halos are almost spherical, if their phase-space DFs 
are approximately isotropic, their DFs could be written as func- 
tions of energy alone: / = f(E). In this case they showed that if 
Q cx r"* 3 , then V(f) would also be described by a power-law, and 
the index a was related to the index f3 through a simple equation. 

In Figure[3]we plot the evolution of V(f) with redshift for all 
the material that lies within twice the virial radius of halo Gi at 2 = 
0. In the bottom right-hand panel, the solid line is a fit to V (/) for 
values of 10~ 4 < / < 10 2 ' 5 . Over this range the best fit power-law 
profile is given by V(f) cx j- 2 - 34 ± 02 which is not very d ifferent 
from the power-law profile fit obtained bv lArad et alj J2004). At all 
redshifts V(f) deviates from the power-law profile at both low and 
high end. It is likely that at the high-/ end the distribution deviates 
from the power-law due to the unresolved subhalos below the mass 
resolution limit of the simulation. 

The Mixing Theorem requires that the coar se-grained vol - 
ume distribution function V(f) always increases jMathuj[l988h . 
We see that V(f) increases steadily from z = 9 onward for 
/ < lM0kpc -3 (kms -1 ) -3 in such a way that as the system 
evolves, there is more volume associated with material with low /. 
At higher values of /, V(f) remains almost constant with decreas- 
ing redshift, probably owing to the fact that this material forms the 
cusps of DM halos. While the Mixing Theorems have previously 
been demonstrated for isolated systems, this is to our knowledge 
the first demonstration of their validity in a cosmological context. 

It is particularly illuminating to plot the full coarse-grained 
phase space density of all particles as a function of radius from 
the center of mass of the main halo at each redshift (Figure gj. At 
any radius from the center there exist particles with phase-space 
densities spanning between 4-8 decades in /. The colored con- 
tours represent a constant particle number per unit area on the plane 
(log(/), log(r)). The yellow/orange contours represent the param- 
eter range with the largest number of particles, while the blue/black 
contours represent regions with the smallest number of particles. 
The solid white curves on each plot show the median value (F) of 
/(x, v) in 100 logarithmically spaced bins in r. These curves are 
identical to the various curves in Figure|2]and trace the regions with 
the largest particle number density at each radius. Numerous spikes 
in / are seen at r > 0.1r v i r and correspond to DM sub-halos. 

As the evolution proceeds, there is an overall lowering of the 



median phase-space density F (white curves). There is also a con- 
tinuous decrease of the low-end envelope of / to lower values 
pointing to an increasingly large amount of matter with low phase- 
space density. However the upper envelope representing the highest 
phase space density regions at the centers of the subhalos remains 
relatively unchanged with redshift, indicating that primordial val- 
ues of / are largely preserved at the centers of DM subhalos lying 
outside the center of the main halo. 

It is illustrative to compare the evolution of F(r) and Q(r) 
with redshift. Figure [5] shows the two curves plotted as a function 
of radius (in comoving units) at each redshift for all the matter in 
halo Gi that lies within two virial radii at z — 0, as a function of 
physical radius from the center of the most massive progenitor of 
the final halo (in units of the virial radius at z — 0). In each panel 
a thin vertical line is drawn at the virial radius of the halo at that 
redshift. 

We find the best fit power-law Q cx r ~ 184±0 012 at z = t o 
the profile within 0.6r v i r - In agreement with Hoffman et aljj2007h . 
we find that the same power-law provides a reasonably good fit to 
the profiles at redshifts from z ~ 5 to z = 0. While Q always 
decreases with radius since it is a spatial average over increas- 
ingly large volumes of configuration space, the median phase-space 
density F decreases monotonically with radius only within about 
0.6r v ir of the main progenitor at that redshift. Outside the virial 
radius at each redshift (i.e. to the right of the vertical dashed line), 
F(r) become significantly flatter, even increasing with increasing 
radius. At all z 7^ this represents the median phase-space density 
of matter that will lie within the outer radius of the halo at z — 0, 
but is either still unvirialized or lies within small halos. In general 
this material is not as mixed as the more massive main progenitor. 
The nearly constant value of F beyond the virial radius at high z 
may be interpreted as the median phase-space density of DM in the 
Universe at that redshift. At lower redshifts, a non-negligible frac- 
tion of this matter lies within inner regions of small halos, which 
have undergone relatively small amount of mixing. As more and 
more material undergoes substantial mixing as the halos evolve, 
the high phase-space density unmixed material in the subhalo cen- 
ters, as well as in the relatively unmixed streams formed from dis- 
rupted subhalos, prevents the median at large radii from decreasing 
rapidly, despite the large increase in low material with low /-values 
within the virialized regions of the main halo. 

To better understand the evolution of / with redshift we plot 
histograms of log(/) at each redshift in Figure [6] (thick solid his- 
tograms). The thin red curves at each redshift represent the best-fit 
Gaussians to the distribution of log(/) at that redshift. A Gaus- 
sian provides a good fit to the majority of the mass. (Since log(/) 
is close to Gaussian, / itself has a log-normal distribution.) By 
approximately z = 1, the skewness of the distribution is signif- 
icant, and the skewness increases steadily until z = 0. To bet- 
ter understand the origin of the matter in the high-/ tail at z=0, 
we plot the distribution of DM particles that have the values of 
/ ^ lA/0kpc~ 3 (kms -1 ) -3 at z = 0. The ordinate values of 
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Figure 3. The volume DF V(f) for Gi halo in the L25 simulation at different redshifts. A power- law fit to V(f) <x z-2.34i0.oi7_ is shown at z = o (solid 
line). 



this distribution (multiplied by a factor of 10 to enhance visibility), 
are shown in the dashed histograms. The thin blue curves show that 
the high-/ tail is also well fitted by a Gaussian (except at z = 0). 
As we saw in Figure [4] most of these high-/ tail particles are in 
subhalos at z — 0, while some are also in the highest phase-space 
density particles in the central cusp of the halo at z = 0. This high- 
/ sub-population has a Gaussian distribution at z = 9 with a mean 
/ = 1.73 ± O.6M0kpc _3 (kms~ 1 )~ 3 compared with the mean 
/ = 1.57 ± O.54M0kpc~ 3 (kms _1 )" 3 for all the particles (in 



the solid curve) at z = 9. A Students' T-test indicated with 99% 
confidence that both distributions are drawn from the same popula- 
tion at z = 9. This implies that the material that lies in the centers 
of DM subhalos at z = will have phase-space densities that are 
representative of the mean phase-space density of matter at z — 9. 
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4 COMPARISON OF FOUR MILKY- WAY SIZED HALOS 

In this section we present results for the phase-space DFs for all 
four of the DM halos described in §[2] Although we present results 
mainly at z — 0, evolution with redshift for each of the halos were 
similar to the evolution of halo Gi presented earlier. 

Figure[7]compares the volume distribution of phase space den- 
sity V(f) for the four different halos at z — 0. All three halos from 



the L25 simulation (Gi, G2, G3) show almost identical profiles in 
V(f) confirming the universality of the process that produced the 
phase-space DF. For halo L20, V(f) lies systematically above the 
other curves especially at higher values of /, where it also extends 
to larg e values of /. It was first shown by ISharma & Stei nmetz 
d2006l) that this is a numerical consequence of the increased mass 
resolution of the simulation. This indicates that the absolute value 
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of / derived in the previous section is somewhat dependent on the 
mass resolution of the simulation and increases slightly with in- 
creasing mass resolution. In all four halos V(f) is well approxi- 
mated by a power-law over nearly six orders of magnitude in /. 
The first column of Table Q] gives details of the power-law fits to 
V(f) (i.e. the power-law slopes and their errors over the range 
10 -4 < f < 10 2 ' 5 .) The power-law sl opes we obtained are similar 
to those obtained by previous authors dArad et alj2004l) . 

In Figure[8]we plot F(r) and Q(r) for the four different halos 
at z = 0. The top set of four curves show Q(r) while the lower 



set of curves show F(r). In all four halos, Q(r) is well fitted by 
a power-law of slope /3 = 1.8 — 1.9 (see Table [TJ. F(r) shows 
significant deviations from a simple power-law profile, with a sys- 
tematic upturn beyond r/r v i r > 0.1. The higher resolution simu- 
lation (L20) appears to have systematically higher F and Q values 
than the halos from the low resolution simula tion (again, a numeri- 
cal co nsequence of the higher mass resolution lSharma & Steinmetj 
l200d) . Tabled gives the values for the slopes and the error-bars on 
the power-law fits to Q(r) and F(r) for r < 0.6rvi r as well as the 
inner power-law slope of F(r). 
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Table 1. Power-law indices for fits to V, Q, and F for four halos at z = for r < 0.6r v j r 



Halo 


V 


Q 


F 


Gi 


-2.34 ±0.02 


-1.84 ±0.01 


-1.59 ±0.05 


G 2 


-2.37 ±0.02 


-1.82 ±0.01 


-1.46 ± 0.06 


G 3 


-2.35 ±0.02 


-1.75 ±0.01 


-1.42 ± 0.03 


L20 


-2.27 ±0.01 


-1.87 ±0.01 


-1.64 ±0.07 
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5 DISCUSSION 

Several previous studies have attempted to account for the ori- 
gin of the power-law Q(r) profiles of DM halos. Notably, it has 
been argued that power-law profiles result from virialization and 
not from the hierarchical sequence of mergers, since they are also 
produced in simple spherical gravitational collapse simulations 



dTavlor & NavarroteOOlkrBarnes et alj200dl2007h . Our results pre- 
sented in the previous section (Figure[5j show that while Q(r) and 
F(r) have approximately power-law form within 0.6r v ir at a given 
redshift, F(r) flattens out and remains quite flat beyond this radius. 
Furthermore, as the hierarchical growth of the halo progresses these 
approximately power-law profiles extend to larger radii until they 
encompass all the mass within the virial radius at z — 0. 

It has been sho wn from both t heoretical argum ents and nu- 
merical simulations (Dehnen 20Q5|; IVass et a l]|2009) that power- 
law profiles of phase-space density for central cusps are well pre- 
served during major mergers. At intermediate times during a major 
merger, deviations from the power-law profile are seen but these 
largely disappear at the end of the merger. In major mergers there 
are two main reasons for the preservation of the power-law pro- 
files. First, the most tightly bound material - that forming a steep 
central cusp or shallow core preserves its phase-space density in the 
final remnant. This is a consequence of the additivity of the excess 
mass function - a resul t of the fact tha t steeper cusps are less mixed 
than shallower cusps dDehnenl l2005). Over 60% of the material 
in the central cusp (within one scale radius) of a pro genitor NFW 
halo r emains within the cusp of the merger remnant dValluri et al.l 
2007). The second reason for preservation of power-law profiles 
at large radii is that material outside three scale radii of the pro- 
genitor halos is redistributed to other radial bins almost uniformly 
from each radial interval. In addition, nearly 40% of the material 
within the virial radius of the progenitor halos is ejected to be- 
yond the virial rad ius of the final remnant dKazantzidis et alj2006l : 
IValluri et al]|2007h . This material can be shown to have originated 
in roughly equal fractions from each radial interval beyond three 
scale radii and has higher phase-space density than expected from 
the simple power-law extrapolation of the inner power-law. This 
self-similar redistribution of material contributes to the preserva- 
tion of power-law profiles in Q(r) and F(r) in equal mass binary 
mergers. Thus, once a power-law phase-space DF has been estab- 
lished in a DM halo, major mergers will not destroy this profile. 
This has be en confirmed by our recent work on mergers of equal 
mass halos dVass et alj|20 09). 

As discussed in the previous section (Figure [6j, the coarse- 
grained phase-space density / in ACDM halos has a nearly log- 
normal distribution with the median of log(/) (the peak of the his- 
togram /peak) evolving steadily toward lower values of / as the 
halo grows. We now quantify the evolution of /peak with redshift. 

Figure [9] shows the evolution of median phase-space density 
/peak derived from the histograms of log(/) in Figure[6]for each 
of the 4 Milky- Way sized cosmological halos in this study. The 
solid dots represent the values of / pea k as a function a, while the 
dotted curves are meant to guide the eye by connecting points for 
each of the 4 individual halos. For a < 0.5 the decline in / pea k 
is approximately power-law: / pca k(ffl) oc a -4 ' 5 . This implies that 
the matter inside the halos has undergone significant mixing due to 
virialization during the hierarchical formation process, the leveling 
off of the curve at a ~ 0.5 indicates a decrease in mixing after 
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Figure 9. The median phase-space density (peak of histograms / pea k m 
Fig. [6) at each redshift as a function of the cosmic scale factor a = (1 + 
z) , plotted for four Milky-way sized halos (points). The thin dotted lines 
connect points for the four individual halos. 



z — 1 , possibly resulting from a decrease in mass accreted in 
major mergers. 

On average, as a halo grows via accretion its median density 
decreases as a power-law with time, despite the fact that the most 
centrally concentrated material retains its original high phase-space 
density. This nearly power-law profile in / pea k leads to some in- 
sights into the development of phase-space density profi les. We 
can b reak down the formation of halos into two phases ihi et al.l 
120071) : the fast accretion regime during which halo mass grows very 
rapidly and the slow accretion regime. In the fast accretion regime, 
the mixing processes are very efficient as the potential well is estab- 
lished and potential fluctuates rapidly and constantly. This results 
in a rapid decrease in the overall central phase-space density of the 
halo (as seen from the rapid drop in the central most regions of 
the profiles in Fig. [2]) The inner profile of phase-space distribution 
should be largely set during this stage. As halos grow subsequently, 
either by major mergers or quiescent accretion of smaller halos, the 
average phase-space density decreases, but the pre-existing high 
central phase-space density cusps are preserved. 

At these later times the evolution of the true phase-space den- 
sity is complex and occurs due to the accretion of high phase-space 
density material in subhalos, as well as loosely bound material at 
the edges of the halo. The steady decrease in the amplitude and 
increase in the slope of the F(r) profiles in Fig. [2] with redshift 
show that the true phase-space density does not obey the power- 
law profiles seen in Q(r) at all redshifts. Variations in the individ- 
ual power-law profiles of halos both at z — (Fig. [8]l and with 



redshift reflect the large variation in the cosmic accretion histories 
of individual halos. Additional deviations from the power-law arise 
due to the matter bound to subhalos that survives with high phase- 
space density and leads to the variation of the coarse-grained / of 
some six orders of magnitude in the outer regions of halos. 

All this indicates that the actual phase-space distribution is 
not as universal and simple as the Q(r) profiles lead one to believe. 
The mechanism behind the universal power-law form of the Q(r) 
profiles is therefore likely to be different (e.g., it is not affected by 
the presence of subhalos) and simpler than the processes that shape 
the distribution of the coarse-grained phase-space density. The Q 
profile is a ratio of the density and velocity dispersion (which can 
be interpreted as a measure of temperature) and is therefore re- 
lated to the entropy. For an ideal monatomic gas the entropy can 
be defined as K gas = T/p 2 ^ 3 , while for DM, by analogy, the en- 
tropy can be defin ed as Kdm = a 1 1 (fJJt dFaltenbacher et al.l2007h . 
This, as noted by iHoffman e t al. (2007), gives Kdm oc Q -2/3 or 
Kdm oc r 1 ' 2 for Q oc r -1 ' 8 . This power-law form and slope of 
the entropy profile is very simila^] to the one found for the gas 
in the outer regions of clusters in co smological simulations (e.g., 
Borg ani et alj|2004 IVoit et al.l 12005). This power-law is also pre- 
dicted by the models of spherical accretion ( Tozzi & Norman! 1200 lh 
and reflects the increasing entropy to which the accreting material 
is heated as the halo grows its mass. Although the processes gov- 
erning the virialization of gas and DM are different (the short-range 
local interactions for the former, and long-scale interactions for the 
latter), the fact that the resulting entropy profiles are quite similar 
indicate that they lead to the same distribution of entropy. The Q(r) 
profile therefore may reflect the overall entropy profile of DM, not 
the coarse-grained local phase-space density, which exhibits a more 
complicated behavior. 



6 SUMMARY AND CONCLUSIONS 

We have investigated the evolution of the phase-space density of the 
DM in cosmological simulations of the formation of Milky- Way- 
sized DM halos. The analysis was carried out using two different 
codes for estimating the phase-space density. Both codes give qual- 
itatively similar results, but the estimated values of phase- space 
density / are quite sensitive to the type of code and, for a given 
code, also depend quite sensitively on the choice of smoothing ker- 
nel used. Based on comparisons of the two codes (FiEstAS and En- 
BiD ) and various smoothing parameters (Appendix), we select the 
EnBiD code with n = 10 kernel smoothing and present results for 
the analysis with this set of parameters. 



2 The slope is even more similar if on e takes into account the r andom bulk 
motions of the gas in estimating K g!ls (Falten bacher et al.l2007t) . 
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The simulations presented in this paper complement our anal- 
ysis o f the evolution o f phase-space density in binary major merg- 
ers dVass et al.l [2009). We confirm that the profiles of Q(r) = 
Pdm/cdm computed by previous authors can be described by a 
power-law Q(r) oc r _1 - 8±0 1 over more than two orders of mag- 
nitude in radius in all halos. The median of the phase-space den- 
sity (F(r)) at given radius r, however, exhibits a more compli- 
cated behavior. Although F(r) is approximately a power-law for 
r < 0.6r v i r , the profiles generally flatten in the outer regions. Sub- 
halos contribute somewhat to this behavior, although their effect 
is limited by the relatively small fraction of mass (< 0.1) bound 
to them. However, in addition to subhalos, a significant fraction of 
high phase-space density matter is in the relatively unmixed ma- 
terial (p_ojisMyinsJreamsJ_^ disrupted subhalos 
(e.g.. lArad etal1l2004l ; Ibiemand et alj|2008l) . The rise in F(r) at 
large radii suggests that the fraction of mass in such streams can be 
substantial. 

This behavior holds at earlier epochs. From z — 5 to z = 0, 
material within r v i T (z) at each redshift follows a power-law in Q 
with an approximate power-law slope of ~ —1.8 to —1.9. In con- 
trast, F(r) can only be well described by a power-law in the in- 
ner regions, and its slope changes continuously with redshift. Be- 
yond the virial radius, Q (a quantity that is obtained by averaging 
over increasingly large volumes) decreases rapidly with radius, but 
the median value of F flattens significantly. We argue that F is a 
more physically meaningful quantity, especially for understanding 
the evolution of phase-space density in collisionless DM halos, as 
it measures the median of the true coarse-grained phase-space den- 
sity. 

At all redshifts, the highest values of phase-space density / are 
found at the centers of DM subhalos. In the center of the main halo, 
the median phase-space density (F) drops by about an order of 
magnitude from z — 9 to z = 0. In contrast, the centers of DM sub- 
halos maintain their high values of / ~ lO 3 M0kpc -3 (kms -1 ) -3 
at all redshifts. The highest values of / at the center of the main 
halo are, therefore, lower and less representative of the primor- 
dial phase-space density of DM particles than the central value 
of / in the high phase-space density subhalos. At r v i r , the de- 
crease in median phase-space density is much more significant, 
with F fa 3OM kpc" 3 (kms -1 )" 3 at z = 9 decreasing to 
F ~ lO _3 M0kpc _3 (kms _1 ) -3 at z — 0, a decrease of over 
four orders of magnitude (Figure[4}. 

The evolution of F(r) and V(f) with redshift are consistent 
with expectations from the Mixing Theorems, which require that 
mixing reduces the overall phase-space density of matter in col- 
lisionless systems and that the volume of phase-space associated 
with any value of / increases due to mixing and relaxation. 

The distribution of / is approximately log-normal until z ~ 3. 
As time progresses, the mean and median of log(/) shift to pro- 
gressively lower values as a larger and larger fraction of matter 
undergoes mixing and moves to lower values of /. Some fraction 
of high phase-space density material does survive in the centers 



of subhalos and in the relatively unmixed streams leftover after 
subhalo disruptions, which skews the distribution. Remarkably, the 
highest phase-space density particles at z = have retained their 
phase-space density since z ~ 9, the earliest epoch we analyzed. 
The phase-space density in the centers of DM subhalos is therefore 
representative of the mean phase-space density of DM at the high- 
est redshifts. This can potentially allow for stronger constraints to 
be plac ed on the nature of DM particles from the Tremain e-Gunn 
bound jTremaine & Gunnlll979l : lHogan & Dalcantonll2000l) . 

The median value of phase-space density decreases with 
decreasing redshift approximately as a power-law described by 
/peak — ffl~ 4,5 . This majority of the decrease in / pca k, is the result 
of mixing within virialized halos which reduces the coarse-grained 
phase-space density of matter that has turned around from the Hub- 
ble flow, much of this mixing occurs prior to z ~ 1, after which the 
rate of mixing in galactic sized halos slows down. 
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APPENDIX A: COMPARISON OF "FiEstAS " AND 
"EnBiD" ANALYSIS OF HALO Gi 

The numerical estimation of coarse-grained phase-space densities 
during the evolution of the four N-body halos presented in this pa- 
per was carried out using tw o publicly availa ble codes "FiEstAS " 
dAscasibar & Binrievl 120051) and "EnBiD " dSharma & Steinmet j 
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l2006b . A compa rison of these two cod es has previously been pre- 
sented by ISharma & Steinmetzl d2006h . who showed for an ana- 
lytic DF (for the spherical Hernquist potential), that "EnBiD with 
n = 10 kernel smoothin g" gave the highes t fidelity to the analyti- 



cal DF. In a related work 



Vass et al] J2009t) confirmed the findings 



of lSharma & Steinmetzl d2006t> for~a spherical isotropic NFW halo 
This latter study is the main basis for our choice of EnBiD with 
n — 10 kernel smoothing. After this p aper was submitted to th e 
Journal we became aware of the work of M acieiewski et al] J2008h . 
These authors carried out a similar, and somewhat more detailed 
comparison, of coarse-grained phase-space density estimators for 
N-body simulations on cosmological halos at z = 0. Our results 
are in agreement with theirs. 

However, it is unclear whether the comparisons with analytic 
profiles of isolated halos at z = are valid for matter distributions 
arising from cosmological N-body distributions, at high redshifts, 
where the majority of the particles actually lie outside virialized 
halos. Our purpose in this Appendix is to present a comparison of 
results obtained with the different codes at a range of redshifts to 
allow readers to appreciate how sensitive some of the results pre- 
sented in this paper are to the choice of code and smoothing pa- 
rameters used in density estimation. Since the coarse-grained DF 
/(x, v) is a 6-dimensional function, it is difficult to compare es- 
timates for this function obtained using different codes. Tradition- 
ally the best single variable function to compare is the volume DF 
V(f). Variations in the estimation of / translate to variations in 
estimation of V(f). In Figure [ATI we compare the volume den- 
sity of phase-space V(f) obtained for the G\ halo at four differ- 
ent redshifts using the FiEstAS code (triple-dot-dashed curves), En- 
BiD code with no smoothing (dot-dashed curves), EnBiD with Fi- 
EstAS smoothing (long-dashed curves), and EnBiD with a n = 10 
kernel (dotted curves). In each case we see that V(f) has a nearly 
power-law distribution over more than six orders of magnitude 



in / at z 



(from 



10" 



-lO^fokpc-^kms -1 )- 3 ). Ta- 



ble IA1I gives the slopes of power-law fits for the four different 
estimates of V(f) over the range 10~ 4 < / < 10 3 at z = 0. 
The FiEstAS estimate of V(f) is systematically lower than all En- 
BiD estimates at high values of / (at all redshifts). In addition all 
EnBiD curves extend to much higher values of / than the FiEs- 
tAS estimate (this is particularly true at z = 9 where the FiEs- 
tAS estimate differs from the other estimates both quantitatively 
and qualitatively). The results from the various EnBiD estimates 
differ very little at intermediate values of / (10 _4 -10 3 ) and conse- 
quently we are confident that conclusions drawn from the median 
/ is quite insensitive to the details of the EnBiD parameters used to 
obtain /. 

The differences between the various estimates of F at large 
radii become significantly larger at higher redshift. In Figure lA2l we 
plot the four different estimates of F(r) at four different redshifts 
in the evolution. The vertical dashed line in each panel represents 
the virial radius of the main halo at that redshift. We see that F from 
any of the codes is quite flat beyond the virial radius in all cases, 
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Figure Al. The volume DF V(f), at different redshifts during the evolu- 
tion of the Gi halo in the L25 simulation. Plots compare results obtained 
with the FiEstAS code and EnBiD code (using three different parameters) 
as indicated in the line legends. 



Table Al. Power-law indices V(f) oc / a and F(r) oc r 
from different codes 



Code 


a 


P 


FiEstAS 


2.62 ±0.06 


1.43 ± 0.02 


EnBiD (no smoothing) 


2.35 ±0.02 


1.65 ± 0.02 


EnBiD (FiEstAS smoothing) 


2.46 ± 0.03 


1.69 ± 0.07 


EnBiD (kernel n = 10) 


2.34 ±0.02 


1.59 ± 0.05 



but the absolute values of the curves differ significantly. Note that 
at z = 9, the different estimates can differ by as much as two orders 
of magnitude. 

In this paper we choose to present results from EnBiD with 
n = 10 kernel smoothing over the other estim ates largely because it 
does the best job of reproducing analytic DFs dSharma & Steinmetzl 
2006) and because it appears to provide a good upper limit to the 
phase-space density both at high and low values of /. In the ab- 
sence of an analytic comparison of the estimates at high redshift, 
we caution the reader to refrain from drawing very strong conclu- 
sions regarding absolute values of / or F from the results presented 
here. 
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